SUBROUTINE CARBON_DIOXIDE &
& (M, Twater, lamcarbdi, ddz, pressure, wind10, h1, l1, ls1, dhw, dhw0, dt, carbdi)

! Subroutine CARBON_DIOXIDE calculates diffusion of carbon dioxide 
! in the water column

use METH_OXYG_CONSTANTS
use PHYS_FUNC, only : &
& HENRY_CONST, &
& GAS_WATATM_FLUX

implicit none

integer(4), intent(in) :: M

real(8), intent(in) :: Twater(1:M+1) ! Water temperature profile, Celsius
real(8), intent(in) :: lamcarbdi(1:M)
real(8), intent(in) :: ddz(1:M)

real(8), intent(in) :: pressure ! Atmospheric pressure, Pa
real(8), intent(in) :: wind10 ! Wind speed at 10 m above the surface
real(8), intent(in) :: h1, l1, ls1
real(8), intent(in) :: dhw, dhw0
real(8), intent(in) :: dt

real(8), intent(inout) :: carbdi(1:M+1,1:2)
integer(4), parameter :: water_carbdi_indic = 10
integer(4), parameter :: vector_length = 350

real(8), allocatable :: a(:), b(:), c(:), f(:), y(:)

real(8), save :: Fcarbdi1 = 0.d0 ! the flux of carbon dioxide at the bottom

real(8) :: x, Flux_atm

allocate (a(1:vector_length),b(1:vector_length),c(1:vector_length), &
&         f(1:vector_length),y(1:vector_length))

if (h1 > 0) then
  if (l1 == 0) then
!    c(1) = 1.
!    b(1) = 0.
!    f(1) = co2_pres_atm0 * &
!    & HENRY_CONST(Henry_const0_co2, Henry_temp_dep_co2, Henry_temp_ref, Twater(1)+273.15) !co2_atm0 
    x = HENRY_CONST(Henry_const0_co2, Henry_temp_dep_co2, Henry_temp_ref, Twater(1)+273.15)
    Flux_atm = GAS_WATATM_FLUX &
    & (Twater(1),wind10,carbdi(1,1),co2_pres_atm0,x,water_carbdi_indic)
  else
    Flux_atm = 0.
  endif
  c(1)   = - lamcarbdi(1)/(h1*ddz(1)) + 0.5*dhw0/dt - & 
  & 0.5*ddz(1)*h1/dt
  b(1)   = - lamcarbdi(1)/(h1*ddz(1)) + 0.5*dhw0/dt
  f(1)   = - 0.5*ddz(1)*h1*carbdi(1,1)/dt + Flux_atm
  call DIFF_COEF(a,b,c,f,2,M,2,M,water_carbdi_indic,dt)
  if (ls1 > 0) then  
    c(M+1) = - lamcarbdi(M)/(ddz(M)*h1)-0.5*ddz(M)*h1/dt + &
    & 0.5*(dhw-dhw0)/dt
    a(M+1) = - lamcarbdi(M)/(ddz(M)*h1)+0.5*(dhw-dhw0)/dt
    f(M+1) = - 0.5*carbdi(M+1,1)*ddz(M)*h1/dt
    call PROGONKA (a,b,c,f,y,1,M+1)
    y(1:M+1) = dmax1(y(1:M+1),0.d0)
    carbdi(1:M+1,2) = y(1:M+1)
  else  
    x = 0.5*ddz(M)*h1/dt
    a(M+1) = - lamcarbdi(M)/(ddz(M)*h1) + 0.5*(dhw-dhw0)/dt
    c(M+1) = a(M+1) - x
    f(M+1) = - carbdi(M+1,1)*x + Fcarbdi1
    call PROGONKA (a,b,c,f,y,1,M+1)
    y(1:M+1) = dmax1(y(1:M+1),0.d0)
    carbdi(1:M+1,2) = y(1:M+1)
  endif
endif 

deallocate (a, b, c, f, y)

END SUBROUTINE CARBON_DIOXIDE